#############################################/
## R-Script Master - If it ain't broke      ##
## by: Christian Schimpf et al.             ##
## -- Replication Material (Analyses File)  ##
## Version: Dec. 16, 2020                   ##
## Last updated: August 25, 2021            ##
## R-Version: 4.1.0.                        ##
#############################################/

###################################################################/
#### /// Section 1: Load data                                 #####
##################################################################/

df_Alberta <- read.csv("./Data_Schimpfetal_Replication_25082021.csv")

###################################################################/
#### /// Section 2: IRT Analysis - Policy Items                #####
##################################################################/

## This section uses a graded response model on the set of the nine
## relevant policy items to score respondents on the their latent
## support for policy positions that facilitate a transition towards 
## less carbon intensive economy. 

## Item Selection for correlation plot:
item_df_set <- c("Q222_1", "Q222_2", "Q222_3", 
                 "Q222_4", "Q222_6", "Q222_7", 
                 "Q222_8", "Q222_9")
item_df <- df_Alberta[item_df_set]
item_df <- as.data.frame(lapply(item_df, as.numeric))

## Get correlation Matrix:
cormat <- round(cor(item_df, use = "complete.obs"),2)
head(cormat)

# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
  cormat[lower.tri(cormat)]<- NA
  return(cormat)
}
upper_tri <- get_upper_tri(cormat)

# reshape data
melted_cormat <- melt(upper_tri, na.rm = TRUE)

# Heatmap
heatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value)) + 
  geom_tile() +
  scale_fill_viridis(discrete=FALSE) +
  theme_light() +
  ylab("") +
  xlab("") +
  theme(legend.position = c(0.2, 0.7),
        axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.25),
        legend.key.height = unit(0.7,"line"),
        legend.key.width = unit(0.3,"line")) + 
  labs(fill = "Correlation")


# Screeplot (Figure B2.2 in the Online Appendix)

#run PCA
pca_obj <- prcomp(drop_na(item_df), scale. = TRUE)

#extract variance explained by each component
var_explained_df <- data.frame(PC= paste0("PC",1:8),
                               var_explained=(pca_obj$sdev)^2/sum((pca_obj$sdev)^2))

#Screeplot (explained variance)
#https://datavizpyr.com/how-to-make-scree-plot-in-r-with-ggplot2/
screeplot <-  ggplot(var_explained_df, aes(x=PC,y=var_explained, group=1))+
  geom_point(size=3)+
  geom_line()+
  theme_light() +
  ylab("Variance Explained") +
  xlab("Principal Components")

## Combine heatmapt and scree plot:
tiff("./Heatmap_Screeplot.tiff",
     units="cm", width=16, height=6.5, res=600, compression = "lzw")
grid.arrange(heatmap, screeplot, nrow = 1)
dev.off()

## Fit GRM:

#fit_irt <- grm(item_df) - Contrained Discrimination Parameters
fit_irt_Constrained <- grm(df_Alberta[,c("Q222_1", "Q222_2", "Q222_3", "Q222_4",
                                         "Q222_6", "Q222_7", "Q222_8", "Q222_9")], constrained = TRUE)

#fit_irt <- grm(item_df) - Unconstrained Discrimination Parameters
fit_irt <- grm(df_Alberta[,c("Q222_1", "Q222_2", "Q222_3", "Q222_4",
                             "Q222_6", "Q222_7", "Q222_8", "Q222_9")])

## Compare constrained vs unconstrained model
anova(fit_irt_Constrained,fit_irt)

# Print summary for text (https://www.r-bloggers.com/2013/01/item-response-modeling-of-customer-satisfaction-the-graded-response-model/)
summary(fit_irt)

## Score respondents and add variable to dataset:
pattern <- ltm::factor.scores(fit_irt, resp.pattern=df_Alberta[,c("Q222_1", "Q222_2", "Q222_3", "Q222_4",
                                                              "Q222_6", "Q222_7", "Q222_8", "Q222_9")])
df_Alberta$IRT_SCORE_Policy <-pattern$score.dat$z1

## Plots: Information Trace Lines
plot(fit_irt, type = "IIC")


## Instrument information and standard error

# First: extract the information values
vals <- plot(fit_irt, type = "IIC", items = 0, plot = FALSE, zrange = c(-6,6))


tiff("./IRT_ScalePlot.tiff",
     units="cm", width=16, height=8, res=600, compression = "lzw")
# Second: generate the plot with the IIC first:
par(mar=c(5, 4, 4, 6) + 0.1)
plot(fit_irt, type = "IIC", item = 0, zrange = c(-6,6), xlim=c(-6, 6),
     xlab = "", ylab = "", main="", col="black")
box()

# Then add the standard error of measurement:
par(new=T)
plot(vals[, "z"], 1 / sqrt(vals[, "test.info"]), type = "l", lwd = 2, lty=2,
     xlab = "Transition latent attitude", ylab = "Information(??)", 
     main = "", xlim=c(-6, 6), axes=F)
# And finally, draw the second y-axis:
mtext("SE(??)",side=4,col="black",line=4) 
axis(4, col="black",col.axis="black",las=0)
legend(2, 19, legend=c("Information(??)", "SE(??)"),
       lty=1:2, cex=0.8)
dev.off()

## Item trace plots

#Extract IIC for each of the 6 items:
l1 <- plot(fit_irt, type = "IIC", items = 1, plot = T, zrange = c(-6,6))
l2 <- plot(fit_irt, type = "IIC", items = 2, plot = T, zrange = c(-6,6))
l3 <- plot(fit_irt, type = "IIC", items = 3, plot = T, zrange = c(-6,6))
l4 <- plot(fit_irt, type = "IIC", items = 4, plot = T, zrange = c(-6,6))
l5 <- plot(fit_irt, type = "IIC", items = 5, plot = T, zrange = c(-6,6))
l6 <- plot(fit_irt, type = "IIC", items = 6, plot = T, zrange = c(-6,6))
l7 <- plot(fit_irt, type = "IIC", items = 7, plot = T, zrange = c(-6,6))
l8 <- plot(fit_irt, type = "IIC", items = 8, plot = T, zrange = c(-6,6))

#convert them into datamframes:
l1df <- as.data.frame(l1)
l2df <- as.data.frame(l2)
l3df <- as.data.frame(l3)
l4df <- as.data.frame(l4)
l5df <- as.data.frame(l5)
l6df <- as.data.frame(l6)
l7df <- as.data.frame(l7)
l8df <- as.data.frame(l8)

# Add all information into one dataframe using l1df as the basis:
l1df$Q222_2 <- l2df$Q222_2
l1df$Q222_3 <- l3df$Q222_3
l1df$Q222_4 <- l4df$Q222_4
l1df$Q222_6 <- l5df$Q222_6
l1df$Q222_7 <- l6df$Q222_7
l1df$Q222_8 <- l7df$Q222_8
l1df$Q222_9 <- l8df$Q222_9

# Reshaping dataset for graphing purposes:
test_data_long <- melt(l1df, id="z")

#Create new variable with ideal name ("Item) for graphing:
test_data_long$Item <- test_data_long$variable

# Assign names to variable "Item":
levels(test_data_long$Item) <- c("Item 1 - CR",
                                 "Item 2 - CT",
                                 "Item 3 - TIEE",
                                 "Item 4 - MFEE",
                                 "Item 5 - MT",
                                 "Item 6 - SSR",
                                 "Item 7 - EOSC",
                                 "Item 8 - BMP")

# Plot the lines:
tiff("./IRT_ItemTracePlot.tiff",
     units="cm", width=16, height=8, res=600, compression = "lzw")
ggplot(data=test_data_long,
       aes(x=z, y=value)) +
  geom_line(aes(linetype=Item, color=Item))+
  theme_light() +
  theme(legend.position=c(0.8,0.75),
        legend.key.height = unit(0.8,"line"),
        legend.key.width = unit(0.8,"line"))+
  xlab("??")+
  ylab("I(??)")+
  ylim(0,3.5)+
  guides(linetype=guide_legend(nrow=4,byrow=TRUE, ""),
         color=guide_legend(nrow=4,byrow=TRUE, ""))+
  scale_color_viridis(discrete=T, begin=0, end=0.75) 
dev.off()



###################################################################/
#### /// Section 3: Descriptive Statistics                    #####
##################################################################/

#### >>> Descriptive statistics of all variables ####

## Generate subset of all the relevant variables:

des_vars <- c("Q222_1", "Q222_2", "Q222_3", "Q222_4", "Q222_6", "Q222_7", "Q222_8", "Q222_9",
              "IRT_SCORE_Policy", "OandG_Fut_Importance", "AB_Dependent_OandG",
   "OandG", "OG_Pride", "ClimateChBelief", "LR", "PID_NDP", "Gender", "Edu", "Age", "Urban")

des_dfAB <- df_Alberta[des_vars]

## Generate descriptives (including NA count)
psych::describe(des_dfAB)

# Table categorical/dummy vars:
table(df_Alberta$Edu)
table(df_Alberta$Urban)
table(df_Alberta$Gender)
table(df_Alberta$ClimateChBelief)
table(df_Alberta$OandG)

#### >>> Descriptive analyses - explanatory variables and outcome variable (IRT Scores) ####

## Independent Variables as Factor Variables:
df_Alberta$OandG_Fut_Importance_Fac <- as.factor(df_Alberta$OandG_Fut_Importance)
levels(df_Alberta$OandG_Fut_Importance_Fac) <- c("Strongly \n disagree", "Disagree", "Agree", "Strongly \n agree")
df_Alberta$AB_Dependent_OandG_Fac <- as.factor(df_Alberta$AB_Dependent_OandG)
levels(df_Alberta$AB_Dependent_OandG_Fac) <- c("Strongly \n disagree", "Disagree", "Agree", "Strongly \n agree")

#Helpful links:
# Add mean to plots: https://datavizpyr.com/add-mean-line-to-ridgeline-plot-in-r-with-ggridges/


# OandG will be most important industry in 25 years
Rainbowplot1 <- 
  df_Alberta %>% 
  filter(!is.na(OandG_Fut_Importance_Fac)) %>%
  ggplot(aes(x = IRT_SCORE_Policy, y = OandG_Fut_Importance_Fac, 
             fill = stat(x))) +
  geom_density_ridges_gradient(scale = 0.9, 
                               rel_min_height = 0.01,
                               alpha = 0.9,
                               quantile_lines=T,
                               quantile_fun=function(x,...)mean(x)) + 
  #low values=strongly oppose to high values=strongly support
  scale_fill_gradient(low = "#440154FF", high = "#FDE725FF") +
  scale_discrete_manual("point_color", values = c("#D55E00"), guide = "none") +
  coord_flip() + 
  theme_light() +
  theme(legend.position = "none",
        axis.title = element_text(size = 12),
        axis.text= element_text(size = 11)) +
  xlab("Support for Energy \n Transition Policies") +
  scale_x_continuous(limits=c(-3,3),
                     breaks = seq(-3,3,1),
                     labels=c("-3" = "-3",
                              "-2" = "-2",
                              "-1" = "-1",
                              "0" = "0",
                              "1" = "1",
                              "2" = "2",
                              "3" = "3"))+
  ylab("'Oil and Gas will be Alberta's most \n important industry in 25 years'")

#Save plot:
tiff("./Rainbowplot1.tiff",
     units="cm", width=16, height=6.5, res=600, compression = "lzw")
Rainbowplot1
dev.off()

# AB too dependent on OandG
Rainbowplot2 <- 
  df_Alberta %>% 
  filter(!is.na(AB_Dependent_OandG_Fac)) %>%
  ggplot(aes(x = IRT_SCORE_Policy, y = AB_Dependent_OandG_Fac, 
             fill = stat(x))) +
  geom_density_ridges_gradient(scale = 0.9, 
                               rel_min_height = 0.01,
                               alpha = 0.9,
                               quantile_lines=T,
                               quantile_fun=function(x,...)mean(x)) + 
  #low values=strongly oppose to high values=strongly support
  scale_fill_gradient(low = "#440154FF", high = "#FDE725FF") +
  coord_flip() +
  theme_light() +
  theme(legend.position = "none",
        axis.title = element_text(size = 12),
        axis.text= element_text(size = 11)) +
  xlab("Support for Energy \n Transition Policies") +
  scale_x_continuous(limits=c(-3,3),
                     breaks = seq(-3,3,1),
                     labels=c("-3" = "-3",
                              "-2" = "-2",
                              "-1" = "-1",
                              "0" = "0",
                              "1" = "1",
                              "2" = "2",
                              "3" = "3"))+
  ylab("'Alberta is too dependent on \n the oil and gas industry'")

#Save plot:
tiff("./Rainbowplot2.tiff",
     units="cm", width=16, height=6.5, res=600, compression = "lzw")
Rainbowplot2
dev.off()

###################################################################/
#### /// Section 4: Regression Analyses - Bayesian Linear Model ####
##################################################################/

#### <<>> Fit the Models ####

# Note on priors: http://mc-stan.org/rstanarm/articles/priors.html (last accessed:
# July 1, 2021)

## Fit bivariate model (no additional covariates)
BAY_model_biv <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance + 
                    AB_Dependent_OandG ,
                    iter = 4000, adapt_delta = 0.99,
                  data=df_Alberta)
## Describe posterior distribution
# (https://easystats.github.io/bayestestR/articles/example1.html)
mcmcReg(BAY_model_biv, format = 'latex', doctype = F, ci=.95, hpdi=T)

## Fit main model (additional covariates included)
BAY_model_main <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance + AB_Dependent_OandG  + OandG + OG_Pride +
                             ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
                           iter = 4000, adapt_delta = 0.99, 
                          data=df_Alberta)
## Describe posterior distribution
mcmcReg(BAY_model_main, doctype = F, ci=.95, hpdi=T)

#### <<>> Assess Convergence ####

# Note: Based on Muth et al 2018.

### Bivariate Model 

# Numerical:
summaryBAY_model_biv <- summary(BAY_model_biv)
print(summaryBAY_model_biv)
# Graphical
plot(BAY_model_biv, "trace", pars = c("OandG_Fut_Importance", "AB_Dependent_OandG"))

# Multivariable Model

# Numerical:
summaryBAY_model_main <- summary(BAY_model_main)
print(summaryBAY_model_main)
# Graphical

tiff("./ConvergencePlot.tiff",
     units="cm", width=16, height=6.5, res=600, compression = "lzw")
plot(BAY_model_main, "trace", pars = c("OandG_Fut_Importance", "AB_Dependent_OandG"))
dev.off()


#### <<>> Assess Model Fit ####

## Bivariate Model
pp_check(BAY_model_biv, nreps = 100)+ xlab("Support for Energy \n Transition Policies")

## Multivariate Model
pp_check(BAY_model_main, nreps = 100)+ xlab("Support for Energy \n Transition Policies")

#### <<>> Plot posterior distributions ####

# Source: http://mjskay.github.io/tidybayes/articles/tidy-rstanarm.html#fit-prediction-curves
#https://mc-stan.org/bayesplot/articles/plotting-mcmc-draws.html
#https://mc-stan.org/bayesplot/

## Plot the posterior distributions of the estimates from the main Model

# Get estimates from 4000 iterations
Est_OandG_Fut_Importance <- BAY_model_main %>%
  spread_draws(`OandG_Fut_Importance`) %>%
  mutate(Estimate = 1)
names(Est_OandG_Fut_Importance)[4] <- "Variable"

Est_AB_Dependent_OandG <- BAY_model_main %>%
  spread_draws(`AB_Dependent_OandG`) %>%
  mutate(Estimate = 2)
names(Est_AB_Dependent_OandG)[4] <- "Variable"

## Combine into dataset:
Estimates_Main <- rbind(Est_OandG_Fut_Importance, Est_AB_Dependent_OandG) 
Estimates_Main$Estimate <- as.factor(Estimates_Main$Estimate)
levels(Estimates_Main$Estimate) <- c("Future Importance of OandG", "AB too dependent on OandG")

# Plot posterior distributions:
tiff("./Posterior_Distributions_MainModel.tiff",
     units="cm", width=16, height=8, res=600, compression = "lzw")
ggplot(Estimates_Main, aes(y = Estimate, x = Variable)) +
  stat_dotsinterval(point_interval= median_hdi, quantiles = 100, .width = c(0.89, 0.95))+
  theme_light()+
  xlab("Effect on support for Energy Transition Policies")
dev.off()


##### /// Robustness Checks 1 - Bayesian implementation - Bivariate Analysis on complete cases ####

# We re-run our bivariate analysis using only those cases that ended up in the final multivariable
# model (Table 2 in the paper). This is to ensure that the differences in estimates between
# the bivariate analysis (Model 1, Table A2) is not a result of the differences in observations
# which resulted from listwise deletion.

#re-run bivariate analysis (for all valid cases)
BAY_model_biv_1 <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance + 
                              AB_Dependent_OandG,
                            iter = 4000, adapt_delta = 0.99,
                            data=df_Alberta)
## Desribe posterior distribution
mcmcReg(BAY_model_biv_1, format = 'latex', doctype = F, ci=.95, hpdi=T)

#re-run bivariate analysis (for all valid cases included in main model after listwise deletion)

#Generate subset of relevant variables
subset_vars <- c("IRT_SCORE_Policy", "OandG_Fut_Importance", "AB_Dependent_OandG", "OandG", "OG_Pride",
            "ClimateChBelief", "LR", "PID_NDP", "Gender", "Edu", "Age", "Urban")
#Generate subset of datafrma
df_Alberta_sub <- df_Alberta[subset_vars]
#Delete cases with missings on any relevant variable
df_Alberta_sub <- df_Alberta_sub[complete.cases(df_Alberta_sub),]

#re-run bivariate analysis
BAY_model_biv_2 <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance + 
                            AB_Dependent_OandG,
                          iter = 4000, adapt_delta = 0.99,
                          data=df_Alberta_sub)
## Desribe posterior distribution
mcmcReg(BAY_model_biv_2, format = 'latex', doctype = F, ci=.95, hpdi=T)


##### /// Robustness Checks 2 - Bayesian implementation - Ordered Logit for all policies ####


#Carbon Regulation (Probit)
fit_bay_CR1 <- brm(
  formula = Q222_1 ~ 1 + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
summary(fit_bay_CR1)
describe_posterior(fit_bay_CR1)
mcmcReg(fit_bay_CR1, doctype = T, ci=.95, hpdi=T)

#Carbon Tax
fit_bay_CT1 <- brm(
  formula = Q222_2 ~ 1 + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
summary(fit_bay_CT1)
describe_posterior(fit_bay_CT1)
mcmcReg(fit_bay_CT1, doctype = T, ci=.95, hpdi=T)

#Tax incentives
fit_bay_TI1 <- brm(
  formula = Q222_3 ~ 1 + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
summary(fit_bay_TI1)
describe_posterior(fit_bay_TI1)
mcmcReg(fit_bay_TI1, doctype = T, ci=.95, hpdi=T)


#Money for energy efficiency
fit_bay_MFEE1 <- brm(
  formula = Q222_4 ~ 1  + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
mcmcReg(fit_bay_MFEE1, doctype = T, ci=.95, hpdi=T)


#Expand mass transit
fit_bay_EMT1 <- brm(
  formula = Q222_6 ~ 1  + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
mcmcReg(fit_bay_EMT1, doctype = T, ci=.95, hpdi=T)


#Subsidize renewables
fit_bay_SRE1 <- brm(
  formula = Q222_7 ~ 1  + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
mcmcReg(fit_bay_SRE1, doctype = T, ci=.95, hpdi=T)

#Increase oil capacity
fit_bay_IOC1 <- brm(
  formula = Q222_8 ~ 1  + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
mcmcReg(fit_bay_IOC1, doctype = T, ci=.95, hpdi=T)


#Build more pipelines
fit_bay_BMP1 <- brm(
  formula = Q222_9 ~ 1  + OandG_Fut_Importance + AB_Dependent_OandG + OandG + OG_Pride +
    ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban,
  data = df_Alberta,
  family = cumulative("probit"),
  inits = "0",
  iter = 2000
)
mcmcReg(fit_bay_BMP1, doctype = T, ci=.95, hpdi=T)


##### /// Robustness Checks 3 - Is the effect driven by the sampling? ####

#### >>> Test 1: Main Variable Effects conditional on education ####

#### <<<>>> Replicate Main Model with interaction between Education and main explanatory variables ####

BAY_RBC_model_Conditional_Education <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance*Edu + AB_Dependent_OandG*Edu + OandG + OG_Pride +
                                      ClimateChBelief + LR + PID_NDP + Gender + Age + Urban,
                                    iter = 4000, adapt_delta = 0.99, 
                                    data=df_Alberta)
mcmcReg(BAY_RBC_model_Conditional_Education, format = 'latex', doctype = F, ci=.95, hpdi=T)


#https://strengejacke.wordpress.com/2017/08/23/going-bayes-rstats/

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_1 <- ggpredict(BAY_RBC_model_Conditional_Education, terms = c("OandG_Fut_Importance", "Edu"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_1, rawdata = TRUE)
# -> no visible difference.

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_2 <- ggpredict(BAY_RBC_model_Conditional_Education, terms = c("AB_Dependent_OandG", "Edu"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_2)
plot(df_RBC3_pred_data_2, rawdata = TRUE)
# -> no visible difference.

#### >>> Test 2: Main Variable Effects conditional on age (in groups) ####

# Generate Age_Group Variable
df_Alberta$Age_Cat <- NA
df_Alberta$Age_Cat[df_Alberta$Age<25] <- 1
df_Alberta$Age_Cat[df_Alberta$Age>24 & df_Alberta$Age<35] <- 2
df_Alberta$Age_Cat[df_Alberta$Age>34 & df_Alberta$Age<45] <- 3
df_Alberta$Age_Cat[df_Alberta$Age>44 & df_Alberta$Age<55] <- 4
df_Alberta$Age_Cat[df_Alberta$Age>54 ] <- 5
df_Alberta$Age_Cat <- as.factor(df_Alberta$Age_Cat)
levels(df_Alberta$Age_Cat) <- c("24 or younger", "25-34", "35-44", "45-54"," 55 or older")

#### <<<>>> Replicate Main Model with interaction between Age_Group and main explanatory variables ####

BAY_RBC_model_Conditional_AgeGroup <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance*Age_Cat + AB_Dependent_OandG*Age_Cat + OandG + OG_Pride +
                                                  ClimateChBelief + LR + PID_NDP + Gender + Edu + Urban,
                                                iter = 4000, adapt_delta = 0.99, 
                                                data=df_Alberta)
mcmcReg(BAY_RBC_model_Conditional_AgeGroup, format = 'latex', doctype = F, ci=.95, hpdi=T)


#https://strengejacke.wordpress.com/2017/08/23/going-bayes-rstats/

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_age_1 <- ggpredict(BAY_RBC_model_Conditional_AgeGroup, terms = c("OandG_Fut_Importance", "Age_Cat"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_age_1)
plot(df_RBC3_pred_data_age_1, rawdata = TRUE)
# -> no visible difference.

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_age_2 <- ggpredict(BAY_RBC_model_Conditional_AgeGroup, terms = c("AB_Dependent_OandG", "Age_Cat"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_age_2)
plot(df_RBC3_pred_data_age_2, rawdata = TRUE)
# -> no visible difference.

#### >>> Test 3: Main Variable Effects conditional on partisanship ####

#### <<<>>> Replicate Main Model with interaction between NDP PiD and main explanatory variables ####

BAY_RBC_model_Conditional_NDP <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance*PID_NDP + AB_Dependent_OandG*PID_NDP + OandG + OG_Pride +
                                                  ClimateChBelief + LR + Gender + Age + Edu + Urban,
                                                iter = 4000, adapt_delta = 0.99, 
                                                data=df_Alberta)
mcmcReg(BAY_RBC_model_Conditional_NDP, format = 'latex', doctype = F, ci=.95, hpdi=T)

#https://strengejacke.wordpress.com/2017/08/23/going-bayes-rstats/

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_NDP_1 <- ggpredict(BAY_RBC_model_Conditional_NDP, terms = c("OandG_Fut_Importance", "PID_NDP"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_NDP_1)
plot(df_RBC3_pred_data_NDP_1, rawdata = TRUE)
# -> some minor difference.

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_NDP_2 <- ggpredict(BAY_RBC_model_Conditional_NDP, terms = c("AB_Dependent_OandG", "PID_NDP"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_NDP_2)
plot(df_RBC3_pred_data_NDP_2, rawdata = TRUE)
# -> no visible difference.

#### >>> Test 4: Main Variable Effects conditional on gender ####

#### <<<>>> Replicate Main Model with interaction between NDP PiD and main explanatory variables ####

BAY_RBC_model_Conditional_Gender <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance*Gender + AB_Dependent_OandG*Gender + OandG + OG_Pride +
                                            ClimateChBelief + LR + PID_NDP + Age + Edu + Urban,
                                          iter = 4000, adapt_delta = 0.99, 
                                          data=df_Alberta)
mcmcReg(BAY_RBC_model_Conditional_Gender, format = 'latex', doctype = F, ci=.95, hpdi=T)


#https://strengejacke.wordpress.com/2017/08/23/going-bayes-rstats/

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_Gender_1 <- ggpredict(BAY_RBC_model_Conditional_Gender, terms = c("OandG_Fut_Importance", "Gender"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_Gender_1)
plot(df_RBC3_pred_data_Gender_1, rawdata = TRUE)
# -> some minor difference.

# Plot interaction effects (predicted values) for "OandG_Fut_Importance*Edu" 
df_RBC3_pred_data_Gender_2 <- ggpredict(BAY_RBC_model_Conditional_Gender, terms = c("AB_Dependent_OandG", "Gender"), ppd=T, mdrt.values = "meansd")
plot(df_RBC3_pred_data_Gender_2)
plot(df_RBC3_pred_data_Gender_2, rawdata = TRUE)
# -> no visible difference.



##### /// Robustness Checks 4 - Additional Covariate (Income) ####

# Note: This result is not shown in the paper or appendix. In the first round
# of comments, a reviewer asked how income might affect the relationship we are
# interested in. In footnote 6 in the paper, we note that income and education
# are highly correlated such that income is likely a proxy for educational effects
# which we adjust for by including the education dummy. The following model 
# corroborates this finding as the results are basically identical to those from
# our primary model (Table 2 in the paper).

BAY_model_main_RBC_Income <- stan_glm(IRT_SCORE_Policy ~ OandG_Fut_Importance + AB_Dependent_OandG  + OandG + OG_Pride +
                             ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban + Income,
                           iter = 4000, adapt_delta = 0.99, 
                           data=df_Alberta)
## Describe posterior distribution
mcmcReg(BAY_model_main_RBC_Income, doctype = F, ci=.95, hpdi=T)

##### /// Alternative Model (Frequentist OLS)  ####

## Re-run primary model using OLS in frequentist approach:

OLS_model_main <- lm(IRT_SCORE_Policy ~ OandG_Fut_Importance + AB_Dependent_OandG  + OandG + OG_Pride +
                             ClimateChBelief + LR + PID_NDP + Gender + Edu + Age + Urban, 
                           data=df_Alberta)
#output for appendix:
summary(OLS_model_main)
################################ Session Info #################################

sessionInfo()

############################### END OF R-Script ################################